function corrMat = correlation_AVAmean(type, seisReal, seisSim, layers, D_mean, Vp_mean, Vs_mean, dSim, vpSim, vsSim, wSeis)
%% COMPUTE THE CORR COEF BETWEEN ALL THE SIMULATIONS WITH THE REAL SEISMIC
%
% type - TYPE OF CORRELATION COMPUTED

% seis_real [cell] - REAL SEISMIC
% seis_sim [cell]  - SIMULATED SEISMIC (just the data from the "noIt" iteration)
% layers           - NUMBER AND SIZE O EACH LAYER USED FOR THE ITERATION
%
% corr_cell [cell] - [size(layes,1) noAngles]
%
% MODIFICATIONS
%   SEP 2019    - major performance update.
%   DEC 2018    - weight for seismic is now input variable and is
%                   automatically calculated for deviations from mean model.
%   DEC 2018    - update on input variable to include iteration number.
%   04 DEC 2015 - MAJOR UPDATE ON PERFORMANCE - got ride of for loops
%   21 APR 2014 - CORRELATION VOLUME AS ONLY size(layers,1) SIZE
%   29 SEP 2013 - BUG FIX. CONSIDER ONLY THE SYNTHETIC AT CURRENT SIMULATION
%   10 JUL 2013 - BUG FIX. WHEN REAL SEISMIC == 0 WITHIN LAYER ind WAS NOT
%                   UPDATING. PERFORMANCE UPDATE
%   05 MAR 2013 - NEW CONDITION TO CHECK IF ORIGINAL SEISMIC IS == 0 IN THE
%               ENTIRE LAYER
%   03 OCT 2012 - LAYERS ARE NOW COMPUTED OUTSIDE IN MAIN
%   11 SEP 2012 - GLOBAL CORRELATION IS COMPUTED CORRECTLY IN A NEW FUNCTION
%   07 SEP 2012 - PEARSON CORRELATION COMPUTED BY HAND
%   04 SEP 2012 - CHANGES FOR SINGLE SIMULATIONS RUN - MEMORY UPDATE.
%   21 AUG 2012 - MEMORY ALLOCATION UPDATE
%   02 AUG 2012 - COMPUTE AND WRITE INTO LOG FILE THE GLOBAL CORRELATION OF
%               SYNTHETIC PER SIMULATION
%   20 JUL 2012 - DEFAULT PARAMETERS FOR MIN AND MAX LAYER SIZE REMOVED;
%               CORRELATION VALUES ARE ALWAYS BTW. 0 AND 1
%   11 JUL 2012 - TYPE ADDED TO RUN WITH QUASI CORRELATION
%   21 JUN 2012 - QUASI CORRELATION ADDED
%   27 APR 2012 - FIXED LAST LAYER IS ALWAYS > MinLayerSize
%   23 APR 2012 - NOW COMPUTES CORR COEF BTW LAYERS WITH VARIABLE SIZE
%   30 MAR 2012 - NOW COMPUTES CORR COEF BTW LAYERS WITH FIXED SIZE
%
% LA - 22 APR 2012
%

%% allocate and initialize
sizeGrid    = size(seisReal);
sizeGather  = size(seisReal{1,1});

seisReal    = cell2mat(seisReal);
seisSim     = cell2mat(seisSim);
corrMat     = zeros(size(layers,1) * sizeGrid(1), size(seisReal,2));
layersAll   = repmat(layers, [sizeGrid(1), 1]);

%% computed the residuals between simulated and mean models (normalized values)
meanDDev    = 1-(abs(dSim-D_mean)   ./ max(dSim(:)));
meanVpDev   = 1-(abs(vpSim-Vp_mean) ./ max(vpSim(:)));
meanVsDev   = 1-(abs(vsSim-Vs_mean) ./ max(vsSim(:)));

meanDev     = (meanVsDev + meanVpDev + meanDDev) ./ 3;
meanDev     = reshape(permute(meanDev, [3 1  2]), [], sizeGrid(2));

clear meanDDev meanVpDev meanVsDev

if type==0 % CORRELATION COEFFICIENT
  
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     Needs revision
%
%
%
%     indLayer    = 1;
%     for n=1:size(layersAll,1)
%         seisRealLayer   = seisReal(indLayer: indLayer+layersAll(n)-1,:);
%         seisSyLayer     = seisSim(indLayer:indLayer+layersAll(n)-1,:);
%         if (sum(seisRealLayer) == 0 & sum(seisSyLayer) == 0)
%             corrMat(n ,:) = 1;
%         else
%             corrMat(n ,:) = (wSeis .* (sum((seisRealLayer - repmat(mean(seisRealLayer),size(seisRealLayer,1),1))...
%                 .*(seisSyLayer - repmat(mean(seisSyLayer),size(seisSyLayer,1),1))))./...
%                 (sqrt(sum((seisRealLayer - repmat(mean(seisRealLayer),size(seisRealLayer,1),1)).^2))...
%                 .*sqrt((sum((seisSyLayer - repmat(mean(seisSyLayer),size(seisSyLayer,1),1)).^2))))) + ((1-wSeis) .*reshape(repmat(mean((meanDev(i,:,indLayer:indLayer+layers(n)-1)),3),size(seis_real{1,1},2),1),1,[]));
%         end
%         indLayer = indLayer+layersAll(n,1);
%     end
    
elseif type==1
    indLayer        = 1;
    for n=1:size(layersAll,1)
        seisRealLayer   = seisReal(indLayer:indLayer+layersAll(n)-1,:);
        seisSyLayer     = seisSim(indLayer:indLayer+layersAll(n)-1,:);
        
        % reorder original matrix and repeat to match seismic on cell2mat
        meanLayer2      = repmat(mean(meanDev(indLayer:indLayer+layersAll(n)-1,:),1,'omitnan'),sizeGather(2),1);
        meanLayer2      = meanLayer2(:)';
        
        corrMat(n,:) = (wSeis .* (2*sum((seisRealLayer) .*seisSyLayer,1,'omitnan')) ...
            ./ (sum(seisRealLayer.^2,1,'omitnan') + sum(seisSyLayer.^2,1,'omitnan')))...
            + ((1 - wSeis) .* meanLayer2);
        
        indLayer = indLayer+layersAll(n,1);
   end   
end

% ensure everything is positive and restore original geometry
corrMat(corrMat < 0) = 0;
corrMat(isnan(corrMat(:,:)))=0;
corrMat = mat2cell(corrMat, repmat(size(layers,1),sizeGrid(1),1),  repmat(sizeGather(2),sizeGrid(2),1));

end
